home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / gamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  60.0 KB  |  1,663 lines
  1. /* specfunc/gamma.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_log.h>
  27. #include <gsl/gsl_sf_psi.h>
  28. #include <gsl/gsl_sf_trig.h>
  29. #include <gsl/gsl_sf_gamma.h>
  30.  
  31. #include "error.h"
  32. #include "check.h"
  33.  
  34. #include "chebyshev.h"
  35. #include "cheb_eval.c"
  36.  
  37. #define LogRootTwoPi_  0.9189385332046727418
  38.  
  39.  
  40. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  41.  
  42. #define FACT_TABLE_MAX  170
  43. #define FACT_TABLE_SIZE (FACT_TABLE_MAX+1)
  44. static struct {int n; double f; long i; } fact_table[FACT_TABLE_SIZE] = {
  45.     { 0,  1.0,     1L     },
  46.     { 1,  1.0,     1L     },
  47.     { 2,  2.0,     2L     },
  48.     { 3,  6.0,     6L     },
  49.     { 4,  24.0,    24L    },
  50.     { 5,  120.0,   120L   },
  51.     { 6,  720.0,   720L   },
  52.     { 7,  5040.0,  5040L  },
  53.     { 8,  40320.0, 40320L },
  54.  
  55.     { 9,  362880.0,     362880L    },
  56.     { 10, 3628800.0,    3628800L   },
  57.     { 11, 39916800.0,   39916800L  },
  58.     { 12, 479001600.0,  479001600L },
  59.  
  60.     { 13, 6227020800.0,                               0 },
  61.     { 14, 87178291200.0,                              0 },
  62.     { 15, 1307674368000.0,                            0 },
  63.     { 16, 20922789888000.0,                           0 },
  64.     { 17, 355687428096000.0,                          0 },
  65.     { 18, 6402373705728000.0,                         0 },
  66.     { 19, 121645100408832000.0,                       0 },
  67.     { 20, 2432902008176640000.0,                      0 },
  68.     { 21, 51090942171709440000.0,                     0 },
  69.     { 22, 1124000727777607680000.0,                   0 },
  70.     { 23, 25852016738884976640000.0,                  0 },
  71.     { 24, 620448401733239439360000.0,                 0 },
  72.     { 25, 15511210043330985984000000.0,               0 },
  73.     { 26, 403291461126605635584000000.0,              0 },
  74.     { 27, 10888869450418352160768000000.0,            0 },
  75.     { 28, 304888344611713860501504000000.0,           0 },
  76.     { 29, 8841761993739701954543616000000.0,          0 },
  77.     { 30, 265252859812191058636308480000000.0,        0 },
  78.     { 31, 8222838654177922817725562880000000.0,       0 },
  79.     { 32, 263130836933693530167218012160000000.0,     0 },
  80.     { 33, 8683317618811886495518194401280000000.0,    0 },
  81.     { 34, 2.95232799039604140847618609644e38,  0 },
  82.     { 35, 1.03331479663861449296666513375e40,  0 },
  83.     { 36, 3.71993326789901217467999448151e41,  0 },
  84.     { 37, 1.37637530912263450463159795816e43,  0 },
  85.     { 38, 5.23022617466601111760007224100e44,  0 },
  86.     { 39, 2.03978820811974433586402817399e46,  0 },
  87.     { 40, 8.15915283247897734345611269600e47,  0 },
  88.     { 41, 3.34525266131638071081700620534e49,  0 },
  89.     { 42, 1.40500611775287989854314260624e51,  0 },
  90.     { 43, 6.04152630633738356373551320685e52,  0 },
  91.     { 44, 2.65827157478844876804362581101e54,  0 },
  92.     { 45, 1.19622220865480194561963161496e56,  0 },
  93.     { 46, 5.50262215981208894985030542880e57,  0 },
  94.     { 47, 2.58623241511168180642964355154e59,  0 },
  95.     { 48, 1.24139155925360726708622890474e61,  0 },
  96.     { 49, 6.08281864034267560872252163321e62,  0 },
  97.     { 50, 3.04140932017133780436126081661e64,  0 },
  98.     { 51, 1.55111875328738228022424301647e66,  0 },
  99.     { 52, 8.06581751709438785716606368564e67,  0 },
  100.     { 53, 4.27488328406002556429801375339e69,  0 },
  101.     { 54, 2.30843697339241380472092742683e71,  0 },
  102.     { 55, 1.26964033536582759259651008476e73,  0 },
  103.     { 56, 7.10998587804863451854045647464e74,  0 },
  104.     { 57, 4.05269195048772167556806019054e76,  0 },
  105.     { 58, 2.35056133128287857182947491052e78,  0 },
  106.     { 59, 1.38683118545689835737939019720e80,  0 },
  107.     { 60, 8.32098711274139014427634118320e81,  0 },
  108.     { 61, 5.07580213877224798800856812177e83,  0 },
  109.     { 62, 3.14699732603879375256531223550e85,  0 },
  110.     { 63, 1.982608315404440064116146708360e87,  0 },
  111.     { 64, 1.268869321858841641034333893350e89,  0 },
  112.     { 65, 8.247650592082470666723170306800e90,  0 },
  113.     { 66, 5.443449390774430640037292402480e92,  0 },
  114.     { 67, 3.647111091818868528824985909660e94,  0 },
  115.     { 68, 2.480035542436830599600990418570e96,  0 },
  116.     { 69, 1.711224524281413113724683388810e98,  0 },
  117.     { 70, 1.197857166996989179607278372170e100,  0 },
  118.     { 71, 8.504785885678623175211676442400e101,  0 },
  119.     { 72, 6.123445837688608686152407038530e103,  0 },
  120.     { 73, 4.470115461512684340891257138130e105,  0 },
  121.     { 74, 3.307885441519386412259530282210e107,  0 },
  122.     { 75, 2.480914081139539809194647711660e109,  0 },
  123.     { 76, 1.885494701666050254987932260860e111,  0 },
  124.     { 77, 1.451830920282858696340707840860e113,  0 },
  125.     { 78, 1.132428117820629783145752115870e115,  0 },
  126.     { 79, 8.946182130782975286851441715400e116,  0 },
  127.     { 80, 7.156945704626380229481153372320e118,  0 },
  128.     { 81, 5.797126020747367985879734231580e120,  0 },
  129.     { 82, 4.753643337012841748421382069890e122,  0 },
  130.     { 83, 3.945523969720658651189747118010e124,  0 },
  131.     { 84, 3.314240134565353266999387579130e126,  0 },
  132.     { 85, 2.817104114380550276949479442260e128,  0 },
  133.     { 86, 2.422709538367273238176552320340e130,  0 },
  134.     { 87, 2.107757298379527717213600518700e132,  0 },
  135.     { 88, 1.854826422573984391147968456460e134,  0 },
  136.     { 89, 1.650795516090846108121691926250e136,  0 },
  137.     { 90, 1.485715964481761497309522733620e138,  0 },
  138.     { 91, 1.352001527678402962551665687590e140,  0 },
  139.     { 92, 1.243841405464130725547532432590e142,  0 },
  140.     { 93, 1.156772507081641574759205162310e144,  0 },
  141.     { 94, 1.087366156656743080273652852570e146,  0 },
  142.     { 95, 1.032997848823905926259970209940e148,  0 },
  143.     { 96, 9.916779348709496892095714015400e149,  0 },
  144.     { 97, 9.619275968248211985332842594960e151,  0 },
  145.     { 98, 9.426890448883247745626185743100e153,  0 },
  146.     { 99, 9.332621544394415268169923885600e155,  0 },
  147.     { 100, 9.33262154439441526816992388563e157,  0 },
  148.     { 101, 9.42594775983835942085162312450e159,  0 },
  149.     { 102, 9.61446671503512660926865558700e161,  0 },
  150.     { 103, 9.90290071648618040754671525458e163,  0 },
  151.     { 104, 1.02990167451456276238485838648e166,  0 },
  152.     { 105, 1.08139675824029090050410130580e168,  0 },
  153.     { 106, 1.146280563734708354534347384148e170,  0 },
  154.     { 107, 1.226520203196137939351751701040e172,  0 },
  155.     { 108, 1.324641819451828974499891837120e174,  0 },
  156.     { 109, 1.443859583202493582204882102460e176,  0 },
  157.     { 110, 1.588245541522742940425370312710e178,  0 },
  158.     { 111, 1.762952551090244663872161047110e180,  0 },
  159.     { 112, 1.974506857221074023536820372760e182,  0 },
  160.     { 113, 2.231192748659813646596607021220e184,  0 },
  161.     { 114, 2.543559733472187557120132004190e186,  0 },
  162.     { 115, 2.925093693493015690688151804820e188,  0 },
  163.     { 116, 3.393108684451898201198256093590e190,  0 },
  164.     { 117, 3.96993716080872089540195962950e192,  0 },
  165.     { 118, 4.68452584975429065657431236281e194,  0 },
  166.     { 119, 5.57458576120760588132343171174e196,  0 },
  167.     { 120, 6.68950291344912705758811805409e198,  0 },
  168.     { 121, 8.09429852527344373968162284545e200,  0 },
  169.     { 122, 9.87504420083360136241157987140e202,  0 },
  170.     { 123, 1.21463043670253296757662432419e205,  0 },
  171.     { 124, 1.50614174151114087979501416199e207,  0 },
  172.     { 125, 1.88267717688892609974376770249e209,  0 },
  173.     { 126, 2.37217324288004688567714730514e211,  0 },
  174.     { 127, 3.01266001845765954480997707753e213,  0 },
  175.     { 128, 3.85620482362580421735677065923e215,  0 },
  176.     { 129, 4.97450422247728744039023415041e217,  0 },
  177.     { 130, 6.46685548922047367250730439554e219,  0 },
  178.     { 131, 8.47158069087882051098456875820e221,  0 },
  179.     { 132, 1.11824865119600430744996307608e224,  0 },
  180.     { 133, 1.48727070609068572890845089118e226,  0 },
  181.     { 134, 1.99294274616151887673732419418e228,  0 },
  182.     { 135, 2.69047270731805048359538766215e230,  0 },
  183.     { 136, 3.65904288195254865768972722052e232,  0 },
  184.     { 137, 5.01288874827499166103492629211e234,  0 },
  185.     { 138, 6.91778647261948849222819828311e236,  0 },
  186.     { 139, 9.61572319694108900419719561353e238,  0 },
  187.     { 140, 1.34620124757175246058760738589e241,  0 },
  188.     { 141, 1.89814375907617096942852641411e243,  0 },
  189.     { 142, 2.69536413788816277658850750804e245,  0 },
  190.     { 143, 3.85437071718007277052156573649e247,  0 },
  191.     { 144, 5.55029383273930478955105466055e249,  0 },
  192.     { 145, 8.04792605747199194484902925780e251,  0 },
  193.     { 146, 1.17499720439091082394795827164e254,  0 },
  194.     { 147, 1.72724589045463891120349865931e256,  0 },
  195.     { 148, 2.55632391787286558858117801578e258,  0 },
  196.     { 149, 3.80892263763056972698595524351e260,  0 },
  197.     { 150, 5.71338395644585459047893286526e262,  0 },
  198.     { 151, 8.62720977423324043162318862650e264,  0 },
  199.     { 152, 1.31133588568345254560672467123e267,  0 },
  200.     { 153, 2.00634390509568239477828874699e269,  0 },
  201.     { 154, 3.08976961384735088795856467036e271,  0 },
  202.     { 155, 4.78914290146339387633577523906e273,  0 },
  203.     { 156, 7.47106292628289444708380937294e275,  0 },
  204.     { 157, 1.17295687942641442819215807155e278,  0 },
  205.     { 158, 1.85327186949373479654360975305e280,  0 },
  206.     { 159, 2.94670227249503832650433950735e282,  0 },
  207.     { 160, 4.71472363599206132240694321176e284,  0 },
  208.     { 161, 7.59070505394721872907517857094e286,  0 },
  209.     { 162, 1.22969421873944943411017892849e289,  0 },
  210.     { 163, 2.00440157654530257759959165344e291,  0 },
  211.     { 164, 3.28721858553429622726333031164e293,  0 },
  212.     { 165, 5.42391066613158877498449501421e295,  0 },
  213.     { 166, 9.00369170577843736647426172359e297,  0 },
  214.     { 167, 1.50361651486499904020120170784e300,  0 },
  215.     { 168, 2.52607574497319838753801886917e302,  0 },
  216.     { 169, 4.26906800900470527493925188890e304,  0 },
  217.     { 170, 7.25741561530799896739672821113e306,  0 },
  218.  
  219.     /*
  220.     { 171, 1.24101807021766782342484052410e309,  0 },
  221.     { 172, 2.13455108077438865629072570146e311,  0 },
  222.     { 173, 3.69277336973969237538295546352e313,  0 },
  223.     { 174, 6.42542566334706473316634250653e315,  0 },
  224.     { 175, 1.12444949108573632830410993864e318,  0 },
  225.     { 176, 1.97903110431089593781523349201e320,  0 },
  226.     { 177, 3.50288505463028580993296328086e322,  0 },
  227.     { 178, 6.23513539724190874168067463993e324,  0 },
  228.     { 179, 1.11608923610630166476084076055e327,  0 },
  229.     { 180, 2.00896062499134299656951336898e329,  0 },
  230.     { 181, 3.63621873123433082379081919786e331,  0 },
  231.     { 182, 6.61791809084648209929929094011e333,  0 },
  232.     { 183, 1.21107901062490622417177024204e336,  0 },
  233.     { 184, 2.22838537954982745247605724535e338,  0 },
  234.     { 185, 4.12251295216718078708070590390e340,  0 },
  235.     { 186, 7.66787409103095626397011298130e342,  0 },
  236.     { 187, 1.43389245502278882136241112750e345,  0 },
  237.     { 188, 2.69571781544284298416133291969e347,  0 },
  238.     { 189, 5.09490667118697324006491921822e349,  0 },
  239.     { 190, 9.68032267525524915612334651460e351,  0 },
  240.     { 191, 1.84894163097375258881955918429e354,  0 },
  241.     { 192, 3.54996793146960497053355363384e356,  0 },
  242.     { 193, 6.85143810773633759312975851330e358,  0 },
  243.     { 194, 1.32917899290084949306717315158e361,  0 },
  244.     { 195, 2.59189903615665651148098764559e363,  0 },
  245.     { 196, 5.08012211086704676250273578535e365,  0 },
  246.     { 197, 1.00078405584080821221303894971e368,  0 },
  247.     { 198, 1.98155243056480026018181712043e370,  0 },
  248.     { 199, 3.94328933682395251776181606966e372,  0 },
  249.     { 200, 7.88657867364790503552363213932e374,  0 }
  250.     */
  251. };
  252.  
  253. #define DOUB_FACT_TABLE_MAX  297
  254. #define DOUB_FACT_TABLE_SIZE (DOUB_FACT_TABLE_MAX+1)
  255. static struct {int n; double f; long i; } doub_fact_table[DOUB_FACT_TABLE_SIZE] = {
  256.   { 0,  1.000000000000000000000000000,    1L    },
  257.   { 1,  1.000000000000000000000000000,    1L    },
  258.   { 2,  2.000000000000000000000000000,    2L    },
  259.   { 3,  3.000000000000000000000000000,    3L    },
  260.   { 4,  8.000000000000000000000000000,    8L    },
  261.   { 5,  15.00000000000000000000000000,    15L   },
  262.   { 6,  48.00000000000000000000000000,    48L    },
  263.   { 7,  105.0000000000000000000000000,    105L  },
  264.   { 8,  384.0000000000000000000000000,    384L  },
  265.   { 9,  945.0000000000000000000000000,    945L  },
  266.   { 10, 3840.000000000000000000000000,    3840L   },
  267.   { 11, 10395.00000000000000000000000,    10395L  },
  268.   { 12, 46080.00000000000000000000000,       46080L       },
  269.   { 13, 135135.0000000000000000000000,       135135L      },
  270.   { 14, 645120.00000000000000000000000,      645120L      },
  271.   { 15, 2.02702500000000000000000000000e6,   2027025L     },
  272.   { 16, 1.03219200000000000000000000000e7,   10321920L    },
  273.   { 17, 3.4459425000000000000000000000e7,    34459425L    },
  274.   { 18, 1.85794560000000000000000000000e8,   185794560L   },
  275.   { 19, 6.5472907500000000000000000000e8,            0 },
  276.   { 20, 3.7158912000000000000000000000e9,            0 },
  277.   { 21, 1.37493105750000000000000000000e10,          0 },
  278.   { 22, 8.1749606400000000000000000000e10,           0 },
  279.   { 23, 3.1623414322500000000000000000e11,           0 },
  280.   { 24, 1.96199055360000000000000000000e12,          0 },
  281.   { 25, 7.9058535806250000000000000000e12,           0 },
  282.   { 26, 5.1011754393600000000000000000e13,           0 },
  283.   { 27, 2.13458046676875000000000000000e14,          0 },
  284.   { 28, 1.42832912302080000000000000000e15,          0 },
  285.   { 29, 6.1902833536293750000000000000e15,           0 },
  286.   { 30, 4.2849873690624000000000000000e16,           0 },
  287.   { 31, 1.91898783962510625000000000000e17,          0 },
  288.   { 32, 1.37119595809996800000000000000e18,          0 },
  289.   { 33, 6.3326598707628506250000000000e18,           0 },
  290.   { 34, 4.6620662575398912000000000000e19,           0 },
  291.   { 35, 2.21643095476699771875000000000e20,          0 },
  292.   { 36, 1.67834385271436083200000000000e21,          0 },
  293.   { 37, 8.2007945326378915593750000000e21,           0 },
  294.   { 38, 6.3777066403145711616000000000e22,           0 },
  295.   { 39, 3.1983098677287777081562500000e23,           0 },
  296.   { 40, 2.55108265612582846464000000000e24,          0 },
  297.   { 41, 1.31130704576879886034406250000e25,          0 },
  298.   { 42, 1.07145471557284795514880000000e26,          0 },
  299.   { 43, 5.6386202968058350994794687500e26,           0 },
  300.   { 44, 4.7144007485205310026547200000e27,           0 },
  301.   { 45, 2.53737913356262579476576093750e28,          0 },
  302.   { 46, 2.16862434431944426122117120000e29,          0 },
  303.   { 47, 1.19256819277443412353990764062e30,          0 },
  304.   { 48, 1.04093968527333324538616217600e31,          0 },
  305.   { 49, 5.8435841445947272053455474391e31,           0 },
  306.   { 50, 5.2046984263666662269308108800e32,           0 },
  307.   { 51, 2.98022791374331087472622919392e33,          0 },
  308.   { 52, 2.70644318171066643800402165760e34,          0 },
  309.   { 53, 1.57952079428395476360490147278e35,          0 },
  310.   { 54, 1.46147931812375987652217169510e36,          0 },
  311.   { 55, 8.6873643685617511998269581003e36,           0 },
  312.   { 56, 8.1842841814930553085241614926e37,           0 },
  313.   { 57, 4.9517976900801981839013661172e38,           0 },
  314.   { 58, 4.7468848252659720789440136657e39,           0 },
  315.   { 59, 2.92156063714731692850180600912e40,       0 },
  316.   { 60, 2.84813089515958324736640819942e41,       0 },
  317.   { 61, 1.78215198865986332638610166557e42,       0 },
  318.   { 62, 1.76584115499894161336717308364e43,       0 },
  319.   { 63, 1.12275575285571389562324404931e44,       0 },
  320.   { 64, 1.13013833919932263255499077353e45,       0 },
  321.   { 65, 7.2979123935621403215510863205e45,        0 },
  322.   { 66, 7.4589130387155293748629391053e46,        0 },
  323.   { 67, 4.8896013036866340154392278347e47,        0 },
  324.   { 68, 5.0720608663265599749067985916e48,        0 },
  325.   { 69, 3.3738248995437774706530672060e49,        0 },
  326.   { 70, 3.5504426064285919824347590141e50,        0 },
  327.   { 71, 2.39541567867608200416367771623e51,       0 },
  328.   { 72, 2.55631867662858622735302649017e52,       0 },
  329.   { 73, 1.74865344543353986303948473285e53,       0 },
  330.   { 74, 1.89167582070515380824123960272e54,       0 },
  331.   { 75, 1.31149008407515489727961354964e55,       0 },
  332.   { 76, 1.43767362373591689426334209807e56,       0 },
  333.   { 77, 1.00984736473786927090530243322e57,       0 },
  334.   { 78, 1.12138542651401517752540683649e58,       0 },
  335.   { 79, 7.9777941814291672401518892225e58,        0 },
  336.   { 80, 8.9710834121121214202032546920e59,        0 },
  337.   { 81, 6.4620132869576254645230302702e60,        0 },
  338.   { 82, 7.3562883979319395645666688474e61,        0 },
  339.   { 83, 5.3634710281748291355541151243e62,        0 },
  340.   { 84, 6.1792822542628292342360018318e63,        0 },
  341.   { 85, 4.5589503739486047652209978556e64,        0 },
  342.   { 86, 5.3141827386660331414429615754e65,        0 },
  343.   { 87, 3.9662868253352861457422681344e66,        0 },
  344.   { 88, 4.6764808100261091644698061863e67,        0 },
  345.   { 89, 3.5299952745484046697106186396e68,        0 },
  346.   { 90, 4.2088327290234982480228255677e69,        0 },
  347.   { 91, 3.2122956998390482494366629620e70,        0 },
  348.   { 92, 3.8721261107016183881809995223e71,        0 },
  349.   { 93, 2.98743500085031487197609655470e72,       0 },
  350.   { 94, 3.6397985440595212848901395509e73,        0 },
  351.   { 95, 2.83806325080779912837729172696e74,       0 },
  352.   { 96, 3.4942066022971404334945339689e75,        0 },
  353.   { 97, 2.75292135328356515452597297515e76,       0 },
  354.   { 98, 3.4243224702511976248246432895e77,        0 },
  355.   { 99, 2.72539213975072950298071324540e78,       0 },
  356.   { 100, 3.4243224702511976248246432895e79,       0 },
  357.   { 101, 2.75264606114823679801052037785e80,      0 },
  358.   { 102, 3.4928089196562215773211361553e81,       0 },
  359.   { 103, 2.83522544298268390195083598919e82,      0 },
  360.   { 104, 3.6325212764424704404139816015e83,       0 },
  361.   { 105, 2.97698671513181809704837778865e84,      0 },
  362.   { 106, 3.8504725530290186668388204976e85,       0 },
  363.   { 107, 3.1853757851910453638417642339e86,       0 },
  364.   { 108, 4.1585103572713401601859261374e87,       0 },
  365.   { 109, 3.4720596058582394465875230149e88,       0 },
  366.   { 110, 4.5743613929984741762045187512e89,       0 },
  367.   { 111, 3.8539861625026457857121505465e90,       0 },
  368.   { 112, 5.1232847601582910773490610013e91,       0 },
  369.   { 113, 4.3550043636279897378547301176e92,       0 },
  370.   { 114, 5.8405446265804518281779295415e93,       0 },
  371.   { 115, 5.0082550181721881985329396352e94,       0 },
  372.   { 116, 6.7750317668333241206863982681e95,       0 },
  373.   { 117, 5.8596583712614601922835393732e96,       0 },
  374.   { 118, 7.9945374848633224624099499564e97,       0 },
  375.   { 119, 6.9729934618011376288174118541e98,       0 },
  376.   { 120, 9.5934449818359869548919399477e99,       0 },
  377.   { 121, 8.4373220887793765308690683435e100,      0 },
  378.   { 122, 1.17040028778399040849681667362e102,       0 },
  379.   { 123, 1.03779061691986331329689540625e103,       0 },
  380.   { 124, 1.45129635685214810653605267528e104,       0 },
  381.   { 125, 1.29723827114982914162111925781e105,       0 },
  382.   { 126, 1.82863340963370661423542637086e106,       0 },
  383.   { 127, 1.64749260436028300985882145742e107,       0 },
  384.   { 128, 2.34065076433114446622134575470e108,       0 },
  385.   { 129, 2.12526545962476508271787968008e109,       0 },
  386.   { 130, 3.04284599363048780608774948111e110,       0 },
  387.   { 131, 2.78409775210844225836042238090e111,       0 },
  388.   { 132, 4.0165567115922439040358293151e112,        0 },
  389.   { 133, 3.7028500103042282036193617666e113,        0 },
  390.   { 134, 5.3821859935336068314080112822e114,        0 },
  391.   { 135, 4.9988475139107080748861383849e115,        0 },
  392.   { 136, 7.3197729512057052907148953438e116,        0 },
  393.   { 137, 6.8484210940576700625940095873e117,        0 },
  394.   { 138, 1.01012866726638733011865555744e119,       0 },
  395.   { 139, 9.5193053207401613870056733264e119,        0 },
  396.   { 140, 1.41418013417294226216611778042e121,       0 },
  397.   { 141, 1.34222205022436275556779993902e122,       0 },
  398.   { 142, 2.00813579052557801227588724819e123,       0 },
  399.   { 143, 1.91937753182083874046195391280e124,       0 },
  400.   { 144, 2.89171553835683233767727763739e125,       0 },
  401.   { 145, 2.78309742114021617366983317355e126,       0 },
  402.   { 146, 4.2219046860009752130088253506e127,        0 },
  403.   { 147, 4.0911532090761177752946547651e128,        0 },
  404.   { 148, 6.2484189352814433152530615189e129,        0 },
  405.   { 149, 6.0958182815234154851890356000e130,        0 },
  406.   { 150, 9.3726284029221649728795922783e131,        0 },
  407.   { 151, 9.2046856051003573826354437561e132,        0 },
  408.   { 152, 1.42463951724416907587769802630e134,       0 },
  409.   { 153, 1.40831689758035467954322289468e135,       0 },
  410.   { 154, 2.19394485655602037685165496051e136,       0 },
  411.   { 155, 2.18289119124954975329199548675e137,       0 },
  412.   { 156, 3.4225539762273917878885817384e138,        0 },
  413.   { 157, 3.4271391702617931126684329142e139,        0 },
  414.   { 158, 5.4076352824392790248639591467e140,        0 },
  415.   { 159, 5.4491512807162510491428083336e141,        0 },
  416.   { 160, 8.6522164519028464397823346347e142,        0 },
  417.   { 161, 8.7731335619531641891199214170e143,        0 },
  418.   { 162, 1.40165906520826112324473821082e145,       0 },
  419.   { 163, 1.43002077059836576282654719098e146,       0 },
  420.   { 164, 2.29872086694154824212137066574e147,       0 },
  421.   { 165, 2.35953427148730350866380286512e148,       0 },
  422.   { 166, 3.8158766391229700819214753051e149,        0 },
  423.   { 167, 3.9404222333837968594685507847e150,        0 },
  424.   { 168, 6.4106727537265897376280785126e151,        0 },
  425.   { 169, 6.6593135744186166925018508262e152,        0 },
  426.   { 170, 1.08981436813352025539677334714e154,       0 },
  427.   { 171, 1.13874262122558345441781649128e155,       0 },
  428.   { 172, 1.87448071318965483928245015709e156,       0 },
  429.   { 173, 1.97002473472025937614282252992e157,       0 },
  430.   { 174, 3.2615964409499994203514632733e158,        0 },
  431.   { 175, 3.4475432857604539082499394274e159,        0 },
  432.   { 176, 5.7404097360719989798185753611e160,        0 },
  433.   { 177, 6.1021516157960034176023927864e161,        0 },
  434.   { 178, 1.02179293302081581840770641427e163,       0 },
  435.   { 179, 1.09228513922748461175082830877e164,       0 },
  436.   { 180, 1.83922727943746847313387154568e165,       0 },
  437.   { 181, 1.97703610200174714726899923887e166,       0 },
  438.   { 182, 3.3473936485761926211036462131e167,        0 },
  439.   { 183, 3.6179760666631972795022686071e168,        0 },
  440.   { 184, 6.1592043133801944228307090322e169,        0 },
  441.   { 185, 6.6932557233269149670791969232e170,        0 },
  442.   { 186, 1.14561200228871616264651187999e172,       0 },
  443.   { 187, 1.25163882026213309884380982464e173,       0 },
  444.   { 188, 2.15375056430278638577544233437e174,       0 },
  445.   { 189, 2.36559737029543155681480056857e175,       0 },
  446.   { 190, 4.0921260721752941329733404353e176,        0 },
  447.   { 191, 4.5182909772642742735162690860e177,        0 },
  448.   { 192, 7.8568820585765647353088136358e178,        0 },
  449.   { 193, 8.7203015861200493478863993359e179,        0 },
  450.   { 194, 1.52423511936385355864990984535e181,       0 },
  451.   { 195, 1.70045880929340962283784787050e182,       0 },
  452.   { 196, 2.98750083395315297495382329688e183,       0 },
  453.   { 197, 3.3499038543080169569905603049e184,        0 },
  454.   { 198, 5.9152516512272428904085701278e185,        0 },
  455.   { 199, 6.6663086700729537444112150067e186,        0 },
  456.   { 200, 1.18305033024544857808171402556e188,       0 },
  457.   { 201, 1.33992804268466370262665421635e189,       0 },
  458.   { 202, 2.38976166709580612772506233164e190,       0 },
  459.   { 203, 2.72005392664986731633210805920e191,       0 },
  460.   { 204, 4.8751138008754445005591271565e192,        0 },
  461.   { 205, 5.5761105496322279984808215214e193,        0 },
  462.   { 206, 1.00427344298034156711518019425e195,       0 },
  463.   { 207, 1.15425488377387119568553005492e196,       0 },
  464.   { 208, 2.08888876139911045959957480403e197,       0 },
  465.   { 209, 2.41239270708739079898275781478e198,       0 },
  466.   { 210, 4.3866663989381319651591070885e199,        0 },
  467.   { 211, 5.0901486119543945858536189892e200,        0 },
  468.   { 212, 9.2997327657488397661373070276e201,        0 },
  469.   { 213, 1.08420165434628604678682084470e203,       0 },
  470.   { 214, 1.99014281187025170995338370390e204,       0 },
  471.   { 215, 2.33103355684451500059166481610e205,       0 },
  472.   { 216, 4.2987084736397436934993088004e206,        0 },
  473.   { 217, 5.0583428183525975512839126509e207,        0 },
  474.   { 218, 9.3711844725346412518284931849e208,        0 },
  475.   { 219, 1.10777707721921886373117687056e210,       0 },
  476.   { 220, 2.06166058395762107540226850068e211,       0 },
  477.   { 221, 2.44818734065447368884590088393e212,       0 },
  478.   { 222, 4.5768864963859187873930360715e213,        0 },
  479.   { 223, 5.4594577696594763261263589712e214,        0 },
  480.   { 224, 1.02522257519044580837604008002e216,       0 },
  481.   { 225, 1.22837799817338217337843076851e217,       0 },
  482.   { 226, 2.31700301993040752692985058084e218,       0 },
  483.   { 227, 2.78841805585357753356903784452e219,       0 },
  484.   { 228, 5.2827668854413291614000593243e220,        0 },
  485.   { 229, 6.3854773479046925518730966640e221,        0 },
  486.   { 230, 1.21503638365150570712201364459e223,       0 },
  487.   { 231, 1.47504526736598397948268532937e224,       0 },
  488.   { 232, 2.81888441007149324052307165546e225,       0 },
  489.   { 233, 3.4368554729627426721946568174e226,        0 },
  490.   { 234, 6.5961895195672941828239876738e227,        0 },
  491.   { 235, 8.0766103614624452796574435210e228,        0 },
  492.   { 236, 1.55670072661788142714646109101e230,       0 },
  493.   { 237, 1.91415665566659953127881411447e231,       0 },
  494.   { 238, 3.7049477293505577966085773966e232,        0 },
  495.   { 239, 4.5748344070431728797563657336e233,        0 },
  496.   { 240, 8.8918745504413387118605857518e234,        0 },
  497.   { 241, 1.10253509209740466402128414180e236,       0 },
  498.   { 242, 2.15183364120680396827026175195e237,       0 },
  499.   { 243, 2.67916027379669333357172046456e238,       0 },
  500.   { 244, 5.2504740845446016825794386748e239,        0 },
  501.   { 245, 6.5639426708018986672507151382e240,        0 },
  502.   { 246, 1.29161662479797201391454191399e242,       0 },
  503.   { 247, 1.62129383968806897081092663913e243,       0 },
  504.   { 248, 3.2032092294989705945080639467e244,        0 },
  505.   { 249, 4.0370216608232917373192073314e245,        0 },
  506.   { 250, 8.0080230737474264862701598667e246,        0 },
  507.   { 251, 1.01329243686664622606712104019e248,       0 },
  508.   { 252, 2.01802181458435147454008028642e249,       0 },
  509.   { 253, 2.56362986527261495194981623168e250,       0 },
  510.   { 254, 5.1257754090442527453318039275e251,        0 },
  511.   { 255, 6.5372561564451681274720313908e252,        0 },
  512.   { 256, 1.31219850471532870280494180544e254,       0 },
  513.   { 257, 1.68007483220640820876031206743e255,       0 },
  514.   { 258, 3.3854721421655480532367498580e256,        0 },
  515.   { 259, 4.3513938154145972606892082546e257,        0 },
  516.   { 260, 8.8022275696304249384155496309e258,        0 },
  517.   { 261, 1.13571378582320988503988335446e260,       0 },
  518.   { 262, 2.30618362324317133386487400329e261,       0 },
  519.   { 263, 2.98692725671504199765489322224e262,       0 },
  520.   { 264, 6.0883247653619723214032673687e263,        0 },
  521.   { 265, 7.9153572302948612937854670389e264,        0 },
  522.   { 266, 1.61949438758628463749326912007e266,       0 },
  523.   { 267, 2.11340038048872796544071969939e267,       0 },
  524.   { 268, 4.3402449587312428284819612418e268,        0 },
  525.   { 269, 5.6850470235146782270355359914e269,        0 },
  526.   { 270, 1.17186613885743556369012953528e271,       0 },
  527.   { 271, 1.54064774337247779952663025366e272,       0 },
  528.   { 272, 3.1874758976922247332371523360e273,        0 },
  529.   { 273, 4.2059683394068643927077005925e274,        0 },
  530.   { 274, 8.7336839596766957690697974006e275,        0 },
  531.   { 275, 1.15664129333688770799461766294e277,       0 },
  532.   { 276, 2.41049677287076803226326408256e278,       0 },
  533.   { 277, 3.2038963825431789511450909263e279,        0 },
  534.   { 278, 6.7011810285807351296918741495e280,        0 },
  535.   { 279, 8.9388709072954692736948036845e281,        0 },
  536.   { 280, 1.87633068800260583631372476186e283,       0 },
  537.   { 281, 2.51182272495002686590823983534e284,       0 },
  538.   { 282, 5.2912525401673484584047038284e285,        0 },
  539.   { 283, 7.1084583116085760305203187340e286,        0 },
  540.   { 284, 1.50271572140752696218693588728e288,       0 },
  541.   { 285, 2.02591061880844416869829083919e289,       0 },
  542.   { 286, 4.2977669632255271118546366376e290,        0 },
  543.   { 287, 5.8143634759802347641640947085e291,        0 },
  544.   { 288, 1.23775688540895180821413535163e293,       0 },
  545.   { 289, 1.68035104455828784684342337075e294,       0 },
  546.   { 290, 3.5894949676859602438209925197e295,        0 },
  547.   { 291, 4.8898215396646176343143620089e296,        0 },
  548.   { 292, 1.04813253056430039119572981576e298,       0 },
  549.   { 293, 1.43271771112173296685410806860e299,       0 },
  550.   { 294, 3.08150963985904315011544565835e300,       0 },
  551.   { 295, 4.2265172478091122522196188024e301,        0 },
  552.   { 296, 9.1212685339827677243417191487e302,        0 },
  553.   { 297, 1.25527562259930633890922678431e304,       0 },
  554.   /*
  555.   { 298, 2.71813802312686478185383230631e305,       0 },
  556.   { 299, 3.7532741115719259533385880851e306,        0 },
  557.   { 300, 8.1544140693805943455614969189e307,  }
  558.   */
  559. };
  560.  
  561.  
  562. /* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1<t<1
  563.  */
  564. static double gstar_a_data[30] = {
  565.   2.16786447866463034423060819465,
  566.  -0.05533249018745584258035832802,
  567.   0.01800392431460719960888319748,
  568.  -0.00580919269468937714480019814,
  569.   0.00186523689488400339978881560,
  570.  -0.00059746524113955531852595159,
  571.   0.00019125169907783353925426722,
  572.  -0.00006124996546944685735909697,
  573.   0.00001963889633130842586440945,
  574.  -6.3067741254637180272515795142e-06,
  575.   2.0288698405861392526872789863e-06,
  576.  -6.5384896660838465981983750582e-07,
  577.   2.1108698058908865476480734911e-07,
  578.  -6.8260714912274941677892994580e-08,
  579.   2.2108560875880560555583978510e-08,
  580.  -7.1710331930255456643627187187e-09,
  581.   2.3290892983985406754602564745e-09,
  582.  -7.5740371598505586754890405359e-10,
  583.   2.4658267222594334398525312084e-10,
  584.  -8.0362243171659883803428749516e-11,
  585.   2.6215616826341594653521346229e-11,
  586.  -8.5596155025948750540420068109e-12,
  587.   2.7970831499487963614315315444e-12,
  588.  -9.1471771211886202805502562414e-13,
  589.   2.9934720198063397094916415927e-13,
  590.  -9.8026575909753445931073620469e-14,
  591.   3.2116773667767153777571410671e-14,
  592.  -1.0518035333878147029650507254e-14,
  593.   3.4144405720185253938994854173e-15,
  594.  -1.0115153943081187052322643819e-15
  595. };
  596. static cheb_series gstar_a_cs = {
  597.   gstar_a_data,
  598.   29,
  599.   -1, 1,
  600.   17
  601. };
  602.  
  603.  
  604. /* Chebyshev coefficients for
  605.  * x^2(Gamma*(x) - 1 - 1/(12x)), x = 4(t+1)+2, -1 < t < 1
  606.  */
  607. static double gstar_b_data[] = {
  608.   0.0057502277273114339831606096782,
  609.   0.0004496689534965685038254147807,
  610.  -0.0001672763153188717308905047405,
  611.   0.0000615137014913154794776670946,
  612.  -0.0000223726551711525016380862195,
  613.   8.0507405356647954540694800545e-06,
  614.  -2.8671077107583395569766746448e-06,
  615.   1.0106727053742747568362254106e-06,
  616.  -3.5265558477595061262310873482e-07,
  617.   1.2179216046419401193247254591e-07,
  618.  -4.1619640180795366971160162267e-08,
  619.   1.4066283500795206892487241294e-08,
  620.  -4.6982570380537099016106141654e-09,
  621.   1.5491248664620612686423108936e-09,
  622.  -5.0340936319394885789686867772e-10,
  623.   1.6084448673736032249959475006e-10,
  624.  -5.0349733196835456497619787559e-11,
  625.   1.5357154939762136997591808461e-11,
  626.  -4.5233809655775649997667176224e-12,
  627.   1.2664429179254447281068538964e-12,
  628.  -3.2648287937449326771785041692e-13,
  629.   7.1528272726086133795579071407e-14,
  630.  -9.4831735252566034505739531258e-15,
  631.  -2.3124001991413207293120906691e-15,
  632.   2.8406613277170391482590129474e-15,
  633.  -1.7245370321618816421281770927e-15,
  634.   8.6507923128671112154695006592e-16,
  635.  -3.9506563665427555895391869919e-16,
  636.   1.6779342132074761078792361165e-16,
  637.  -6.0483153034414765129837716260e-17
  638. };
  639. static cheb_series gstar_b_cs = {
  640.   gstar_b_data,
  641.   29,
  642.   -1, 1,
  643.   18
  644. };
  645.  
  646.  
  647. /* coefficients for gamma=7, kmax=8  Lanczos method */
  648. static double lanczos_7_c[9] = {
  649.   0.99999999999980993227684700473478,
  650.   676.520368121885098567009190444019,
  651.  -1259.13921672240287047156078755283,
  652.   771.3234287776530788486528258894,
  653.  -176.61502916214059906584551354,
  654.   12.507343278686904814458936853,
  655.  -0.13857109526572011689554707,
  656.   9.984369578019570859563e-6,
  657.   1.50563273514931155834e-7
  658. };
  659.  
  660. /* complex version of Lanczos method; this is not safe for export
  661.  * since it becomes bad in the left half-plane
  662.  */
  663. static
  664. int
  665. lngamma_lanczos_complex(double zr, double zi, gsl_sf_result * yr, gsl_sf_result * yi)
  666. {
  667.   int k;
  668.   gsl_sf_result log1_r,    log1_i;
  669.   gsl_sf_result logAg_r,   logAg_i;
  670.   double Ag_r, Ag_i;
  671.   double yi_tmp_val, yi_tmp_err;
  672.  
  673.   zr -= 1.0; /* Lanczos writes z! instead of Gamma(z) */
  674.  
  675.   Ag_r = lanczos_7_c[0];
  676.   Ag_i = 0.0;
  677.   for(k=1; k<=8; k++) {
  678.     double R = zr + k;
  679.     double I = zi;
  680.     double a = lanczos_7_c[k] / (R*R + I*I);
  681.     Ag_r +=  a * R;
  682.     Ag_i -=  a * I;
  683.   }
  684.  
  685.   gsl_sf_complex_log_e(zr + 7.5, zi, &log1_r,  &log1_i);
  686.   gsl_sf_complex_log_e(Ag_r, Ag_i,   &logAg_r, &logAg_i);
  687.  
  688.   /* (z+0.5)*log(z+7.5) - (z+7.5) + LogRootTwoPi_ + log(Ag(z)) */
  689.   yr->val = (zr+0.5)*log1_r.val - zi*log1_i.val - (zr+7.5) + LogRootTwoPi_ + logAg_r.val;
  690.   yi->val = zi*log1_r.val + (zr+0.5)*log1_i.val - zi + logAg_i.val;
  691.   yr->err = 4.0 * GSL_DBL_EPSILON * fabs(yr->val);
  692.   yi->err = 4.0 * GSL_DBL_EPSILON * fabs(yi->val);
  693.   yi_tmp_val = yi->val;
  694.   yi_tmp_err = yi->err;
  695.   gsl_sf_angle_restrict_symm_err_e(yi_tmp_val, yi);
  696.   yi->err += yi_tmp_err;
  697.   return GSL_SUCCESS;
  698. }
  699.  
  700.  
  701. /* Lanczos method for real x > 0;
  702.  * gamma=7, truncated at 1/(z+8) 
  703.  * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86]
  704.  */
  705. static
  706. int
  707. lngamma_lanczos(double x, gsl_sf_result * result)
  708. {
  709.   int k;
  710.   double Ag;
  711.   double term1, term2;
  712.  
  713.   x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */
  714.  
  715.   Ag = lanczos_7_c[0];
  716.   for(k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); }
  717.  
  718.   /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */
  719.   term1 = (x+0.5)*log((x+7.5)/M_E);
  720.   term2 = LogRootTwoPi_ + log(Ag);
  721.   result->val  = term1 + (term2 - 7.0);
  722.   result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + 7.0);
  723.   result->err += GSL_DBL_EPSILON * fabs(result->val);
  724.  
  725.   return GSL_SUCCESS;
  726. }
  727.  
  728.  
  729. /* x = eps near zero
  730.  * gives double-precision for |eps| < 0.02
  731.  */
  732. static
  733. int
  734. lngamma_sgn_0(double eps, gsl_sf_result * lng, double * sgn)
  735. {
  736.   /* calculate series for g(eps) = Gamma(eps) eps - 1/(1+eps) - eps/2 */
  737.   const double c1  = -0.07721566490153286061;
  738.   const double c2  = -0.01094400467202744461;
  739.   const double c3  =  0.09252092391911371098;
  740.   const double c4  = -0.01827191316559981266;
  741.   const double c5  =  0.01800493109685479790;
  742.   const double c6  = -0.00685088537872380685;
  743.   const double c7  =  0.00399823955756846603;
  744.   const double c8  = -0.00189430621687107802;
  745.   const double c9  =  0.00097473237804513221;
  746.   const double c10 = -0.00048434392722255893;
  747.   const double g6  = c6+eps*(c7+eps*(c8 + eps*(c9 + eps*c10)));
  748.   const double g   = eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*g6)))));
  749.  
  750.   /* calculate Gamma(eps) eps, a positive quantity */
  751.   const double gee = g + 1.0/(1.0+eps) + 0.5*eps;
  752.  
  753.   lng->val = log(gee/fabs(eps));
  754.   lng->err = 4.0 * GSL_DBL_EPSILON * fabs(lng->val);
  755.   *sgn = GSL_SIGN(eps);
  756.  
  757.   return GSL_SUCCESS;
  758. }
  759.  
  760.  
  761. /* x near a negative integer
  762.  * Calculates sign as well as log(|gamma(x)|).
  763.  * x = -N + eps
  764.  * assumes N >= 1
  765.  */
  766. static
  767. int
  768. lngamma_sgn_sing(int N, double eps, gsl_sf_result * lng, double * sgn)
  769. {
  770.   if(eps == 0.0) {
  771.     lng->val = 0.0;
  772.     lng->err = 0.0;
  773.     *sgn = 0.0;
  774.     GSL_ERROR ("error", GSL_EDOM);
  775.   }
  776.   else if(N == 1) {
  777.     /* calculate series for
  778.      * g = eps gamma(-1+eps) + 1 + eps/2 (1+3eps)/(1-eps^2)
  779.      * double-precision for |eps| < 0.02
  780.      */
  781.     const double c0 =  0.07721566490153286061;
  782.     const double c1 =  0.08815966957356030521;
  783.     const double c2 = -0.00436125434555340577;
  784.     const double c3 =  0.01391065882004640689;
  785.     const double c4 = -0.00409427227680839100;
  786.     const double c5 =  0.00275661310191541584;
  787.     const double c6 = -0.00124162645565305019;
  788.     const double c7 =  0.00065267976121802783;
  789.     const double c8 = -0.00032205261682710437;
  790.     const double c9 =  0.00016229131039545456;
  791.     const double g5 = c5 + eps*(c6 + eps*(c7 + eps*(c8 + eps*c9)));
  792.     const double g  = eps*(c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*g5)))));
  793.  
  794.     /* calculate eps gamma(-1+eps), a negative quantity */
  795.     const double gam_e = g - 1.0 - 0.5*eps*(1.0+3.0*eps)/(1.0 - eps*eps);
  796.  
  797.     lng->val = log(fabs(gam_e)/fabs(eps));
  798.     lng->err = 2.0 * GSL_DBL_EPSILON * fabs(lng->val);
  799.     *sgn = ( eps > 0.0 ? -1.0 : 1.0 );
  800.     return GSL_SUCCESS;
  801.   }
  802.   else {
  803.     double g;
  804.  
  805.     /* series for sin(Pi(N+1-eps))/(Pi eps) modulo the sign
  806.      * double-precision for |eps| < 0.02
  807.      */
  808.     const double cs1 = -1.6449340668482264365;
  809.     const double cs2 =  0.8117424252833536436;
  810.     const double cs3 = -0.1907518241220842137;
  811.     const double cs4 =  0.0261478478176548005;
  812.     const double cs5 = -0.0023460810354558236;
  813.     const double e2  = eps*eps;
  814.     const double sin_ser = 1.0 + e2*(cs1+e2*(cs2+e2*(cs3+e2*(cs4+e2*cs5))));
  815.  
  816.     /* calculate series for ln(gamma(1+N-eps))
  817.      * double-precision for |eps| < 0.02
  818.      */
  819.     double aeps = fabs(eps);
  820.     double c1, c2, c3, c4, c5, c6, c7;
  821.     double lng_ser;
  822.     gsl_sf_result c0;
  823.     gsl_sf_result psi_0;
  824.     gsl_sf_result psi_1;
  825.     gsl_sf_result psi_2;
  826.     gsl_sf_result psi_3;
  827.     gsl_sf_result psi_4;
  828.     gsl_sf_result psi_5;
  829.     gsl_sf_result psi_6;
  830.     psi_2.val = 0.0;
  831.     psi_3.val = 0.0;
  832.     psi_4.val = 0.0;
  833.     psi_5.val = 0.0;
  834.     psi_6.val = 0.0;
  835.     gsl_sf_lnfact_e(N, &c0);
  836.     gsl_sf_psi_int_e(N+1, &psi_0);
  837.     gsl_sf_psi_1_int_e(N+1, &psi_1);
  838.     if(aeps > 0.00001) gsl_sf_psi_n_e(2, N+1.0, &psi_2);
  839.     if(aeps > 0.0002)  gsl_sf_psi_n_e(3, N+1.0, &psi_3);
  840.     if(aeps > 0.001)   gsl_sf_psi_n_e(4, N+1.0, &psi_4);
  841.     if(aeps > 0.005)   gsl_sf_psi_n_e(5, N+1.0, &psi_5);
  842.     if(aeps > 0.01)    gsl_sf_psi_n_e(6, N+1.0, &psi_6);
  843.     c1 = psi_0.val;
  844.     c2 = psi_1.val/2.0;
  845.     c3 = psi_2.val/6.0;
  846.     c4 = psi_3.val/24.0;
  847.     c5 = psi_4.val/120.0;
  848.     c6 = psi_5.val/720.0;
  849.     c7 = psi_6.val/5040.0;
  850.     lng_ser = c0.val-eps*(c1-eps*(c2-eps*(c3-eps*(c4-eps*(c5-eps*(c6-eps*c7))))));
  851.  
  852.     /* calculate
  853.      * g = ln(|eps gamma(-N+eps)|)
  854.      *   = -ln(gamma(1+N-eps)) + ln(|eps Pi/sin(Pi(N+1+eps))|)
  855.      */
  856.     g = -lng_ser - log(sin_ser);
  857.  
  858.     lng->val = g - log(fabs(eps));
  859.     lng->err = c0.err + 2.0 * GSL_DBL_EPSILON * (fabs(g) + fabs(lng->val));
  860.  
  861.     *sgn = ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * ( eps > 0.0 ? 1.0 : -1.0 );
  862.  
  863.     return GSL_SUCCESS;
  864.   }
  865. }
  866.  
  867.  
  868. /* This gets bad near the negative half axis. However, this
  869.  * region can be avoided by use of the reflection formula, as usual.
  870.  * Only the first two terms of the series are kept.
  871.  */
  872. #if 0
  873. static
  874. int
  875. lngamma_complex_stirling(const double zr, const double zi, double * lg_r, double * arg)
  876. {
  877.   double re_zinv,  im_zinv;
  878.   double re_zinv2, im_zinv2;
  879.   double re_zinv3, im_zinv3;
  880.   double re_zhlnz, im_zhlnz;
  881.   double r, lnr, theta;
  882.   gsl_sf_complex_log_e(zr, zi, &lnr, &theta);  /* z = r e^{i theta} */
  883.   r = exp(lnr);
  884.   re_zinv =  (zr/r)/r;
  885.   im_zinv = -(zi/r)/r;
  886.   re_zinv2 = re_zinv*re_zinv - im_zinv*im_zinv;
  887.   re_zinv2 = 2.0*re_zinv*im_zinv;
  888.   re_zinv3 = re_zinv2*re_zinv - im_zinv2*im_zinv;
  889.   re_zinv3 = re_zinv2*im_zinv + im_zinv2*re_zinv;
  890.   re_zhlnz = (zr - 0.5)*lnr - zi*theta;
  891.   im_zhlnz = zi*lnr + zr*theta;
  892.   *lg_r = re_zhlnz - zr + 0.5*(M_LN2+M_LNPI) + re_zinv/12.0 - re_zinv3/360.0;
  893.   *arg  = im_zhlnz - zi + 1.0/12.0*im_zinv - im_zinv3/360.0;
  894.   return GSL_SUCCESS;
  895. }
  896. #endif /* 0 */
  897.  
  898.  
  899. inline
  900. static
  901. int
  902. lngamma_1_pade(const double eps, gsl_sf_result * result)
  903. {
  904.   /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps
  905.    * plus a correction series.
  906.    */
  907.   const double n1 = -1.0017419282349508699871138440;
  908.   const double n2 =  1.7364839209922879823280541733;
  909.   const double d1 =  1.2433006018858751556055436011;
  910.   const double d2 =  5.0456274100274010152489597514;
  911.   const double num = (eps + n1) * (eps + n2);
  912.   const double den = (eps + d1) * (eps + d2);
  913.   const double pade = 2.0816265188662692474880210318 * num / den;
  914.   const double c0 =  0.004785324257581753;
  915.   const double c1 = -0.01192457083645441;
  916.   const double c2 =  0.01931961413960498;
  917.   const double c3 = -0.02594027398725020;
  918.   const double c4 =  0.03141928755021455;
  919.   const double eps5 = eps*eps*eps*eps*eps;
  920.   const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps))));
  921.   result->val = eps * (pade + corr);
  922.   result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  923.   return GSL_SUCCESS;
  924. }
  925.  
  926. inline
  927. static
  928. int
  929. lngamma_2_pade(const double eps, gsl_sf_result * result)
  930. {
  931.   /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps
  932.    * plus a correction series.
  933.    */
  934.   const double n1 = 1.000895834786669227164446568;
  935.   const double n2 = 4.209376735287755081642901277;
  936.   const double d1 = 2.618851904903217274682578255;
  937.   const double d2 = 10.85766559900983515322922936;
  938.   const double num = (eps + n1) * (eps + n2);
  939.   const double den = (eps + d1) * (eps + d2);
  940.   const double pade = 2.85337998765781918463568869 * num/den;
  941.   const double c0 =  0.0001139406357036744;
  942.   const double c1 = -0.0001365435269792533;
  943.   const double c2 =  0.0001067287169183665;
  944.   const double c3 = -0.0000693271800931282;
  945.   const double c4 =  0.0000407220927867950;
  946.   const double eps5 = eps*eps*eps*eps*eps;
  947.   const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps))));
  948.   result->val = eps * (pade + corr);
  949.   result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  950.   return GSL_SUCCESS;
  951. }
  952.  
  953.  
  954. /* series for gammastar(x)
  955.  * double-precision for x > 10.0
  956.  */
  957. static
  958. int
  959. gammastar_ser(const double x, gsl_sf_result * result)
  960. {
  961.   /* Use the Stirling series for the correction to Log(Gamma(x)),
  962.    * which is better behaved and easier to compute than the
  963.    * regular Stirling series for Gamma(x). 
  964.    */
  965.   const double y = 1.0/(x*x);
  966.   const double c0 =  1.0/12.0;
  967.   const double c1 = -1.0/360.0;
  968.   const double c2 =  1.0/1260.0;
  969.   const double c3 = -1.0/1680.0;
  970.   const double c4 =  1.0/1188.0;
  971.   const double c5 = -691.0/360360.0;
  972.   const double c6 =  1.0/156.0;
  973.   const double c7 = -3617.0/122400.0;
  974.   const double ser = c0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7))))));
  975.   result->val = exp(ser/x);
  976.   result->err = 2.0 * GSL_DBL_EPSILON * result->val * GSL_MAX_DBL(1.0, ser/x);
  977.   return GSL_SUCCESS;
  978. }
  979.  
  980.  
  981. /* Chebyshev expansion for log(gamma(x)/gamma(8))
  982.  * 5 < x < 10
  983.  * -1 < t < 1
  984.  */
  985. static double gamma_5_10_data[24] = {
  986.  -1.5285594096661578881275075214,
  987.   4.8259152300595906319768555035,
  988.   0.2277712320977614992970601978,
  989.  -0.0138867665685617873604917300,
  990.   0.0012704876495201082588139723,
  991.  -0.0001393841240254993658962470,
  992.   0.0000169709242992322702260663,
  993.  -2.2108528820210580075775889168e-06,
  994.   3.0196602854202309805163918716e-07,
  995.  -4.2705675000079118380587357358e-08,
  996.   6.2026423818051402794663551945e-09,
  997.  -9.1993973208880910416311405656e-10,
  998.   1.3875551258028145778301211638e-10,
  999.  -2.1218861491906788718519522978e-11,
  1000.   3.2821736040381439555133562600e-12,
  1001.  -5.1260001009953791220611135264e-13,
  1002.   8.0713532554874636696982146610e-14,
  1003.  -1.2798522376569209083811628061e-14,
  1004.   2.0417711600852502310258808643e-15,
  1005.  -3.2745239502992355776882614137e-16,
  1006.   5.2759418422036579482120897453e-17,
  1007.  -8.5354147151695233960425725513e-18,
  1008.   1.3858639703888078291599886143e-18,
  1009.  -2.2574398807738626571560124396e-19
  1010. };
  1011. static const cheb_series gamma_5_10_cs = {
  1012.   gamma_5_10_data,
  1013.   23,
  1014.   -1, 1,
  1015.   11
  1016. };
  1017.  
  1018.  
  1019. /* gamma(x) for x >= 1/2
  1020.  * assumes x >= 1/2
  1021.  */
  1022. static
  1023. int
  1024. gamma_xgthalf(const double x, gsl_sf_result * result)
  1025. {
  1026.   /* CHECK_POINTER(result) */
  1027.  
  1028.   if(x == 0.5) {
  1029.     result->val = 1.77245385090551602729817;
  1030.     result->err = GSL_DBL_EPSILON * result->val;
  1031.     return GSL_SUCCESS;
  1032.   } else if (x <= (FACT_TABLE_MAX + 1.0) && x == floor(x)) {
  1033.     int n = (int) floor (x);
  1034.     result->val = fact_table[n - 1].f;
  1035.     result->err = GSL_DBL_EPSILON * result->val;
  1036.     return GSL_SUCCESS;
  1037.   }    
  1038.   else if(fabs(x - 1.0) < 0.01) {
  1039.     /* Use series for Gamma[1+eps] - 1/(1+eps).
  1040.      */
  1041.     const double eps = x - 1.0;
  1042.     const double c1 =  0.4227843350984671394;
  1043.     const double c2 = -0.01094400467202744461;
  1044.     const double c3 =  0.09252092391911371098;
  1045.     const double c4 = -0.018271913165599812664;
  1046.     const double c5 =  0.018004931096854797895;
  1047.     const double c6 = -0.006850885378723806846;
  1048.     const double c7 =  0.003998239557568466030;
  1049.     result->val = 1.0/x + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*c7))))));
  1050.     result->err = GSL_DBL_EPSILON;
  1051.     return GSL_SUCCESS;
  1052.   }
  1053.   else if(fabs(x - 2.0) < 0.01) {
  1054.     /* Use series for Gamma[1 + eps].
  1055.      */
  1056.     const double eps = x - 2.0;
  1057.     const double c1 =  0.4227843350984671394;
  1058.     const double c2 =  0.4118403304264396948;
  1059.     const double c3 =  0.08157691924708626638;
  1060.     const double c4 =  0.07424901075351389832;
  1061.     const double c5 = -0.00026698206874501476832;
  1062.     const double c6 =  0.011154045718130991049;
  1063.     const double c7 = -0.002852645821155340816;
  1064.     const double c8 =  0.0021039333406973880085;
  1065.     result->val = 1.0 + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8)))))));
  1066.     result->err = GSL_DBL_EPSILON;
  1067.     return GSL_SUCCESS;
  1068.   }
  1069.   else if(x < 5.0) {
  1070.     /* Exponentiating the logarithm is fine, as
  1071.      * long as the exponential is not so large
  1072.      * that it greatly amplifies the error.
  1073.      */
  1074.     gsl_sf_result lg;
  1075.     lngamma_lanczos(x, &lg);
  1076.     result->val = exp(lg.val);
  1077.     result->err = result->val * (lg.err + 2.0 * GSL_DBL_EPSILON);
  1078.     return GSL_SUCCESS;
  1079.   }
  1080.   else if(x < 10.0) {
  1081.     /* This is a sticky area. The logarithm
  1082.      * is too large and the gammastar series
  1083.      * is not good.
  1084.      */
  1085.     const double gamma_8 = 5040.0;
  1086.     const double t = (2.0*x - 15.0)/5.0;
  1087.     gsl_sf_result c;
  1088.     cheb_eval_e(&gamma_5_10_cs, t, &c);
  1089.     result->val  = exp(c.val) * gamma_8;
  1090.     result->err  = result->val * c.err;
  1091.     result->err += 2.0 * GSL_DBL_EPSILON * result->val;
  1092.     return GSL_SUCCESS;
  1093.   }
  1094.   else if(x < GSL_SF_GAMMA_XMAX) {
  1095.     /* We do not want to exponentiate the logarithm
  1096.      * if x is large because of the inevitable
  1097.      * inflation of the error. So we carefully
  1098.      * use pow() and exp() with exact quantities.
  1099.      */
  1100.     double p = pow(x, 0.5*x);
  1101.     double e = exp(-x);
  1102.     double q = (p * e) * p;
  1103.     double pre = M_SQRT2 * M_SQRTPI * q/sqrt(x);
  1104.     gsl_sf_result gstar;
  1105.     int stat_gs = gammastar_ser(x, &gstar);
  1106.     result->val = pre * gstar.val;
  1107.     result->err = (x + 2.5) * GSL_DBL_EPSILON * result->val;
  1108.     return stat_gs;
  1109.   }
  1110.   else {
  1111.     OVERFLOW_ERROR(result);
  1112.   }
  1113. }
  1114.  
  1115.  
  1116. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  1117.  
  1118.  
  1119. int gsl_sf_lngamma_e(double x, gsl_sf_result * result)
  1120. {
  1121.   /* CHECK_POINTER(result) */
  1122.  
  1123.   if(fabs(x - 1.0) < 0.01) {
  1124.     /* Note that we must amplify the errors
  1125.      * from the Pade evaluations because of
  1126.      * the way we must pass the argument, i.e.
  1127.      * writing (1-x) is a loss of precision
  1128.      * when x is near 1.
  1129.      */
  1130.     int stat = lngamma_1_pade(x - 1.0, result);
  1131.     result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 1.0));
  1132.     return stat;
  1133.   }
  1134.   else if(fabs(x - 2.0) < 0.01) {
  1135.     int stat = lngamma_2_pade(x - 2.0, result);
  1136.     result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 2.0));
  1137.     return stat;
  1138.   }
  1139.   else if(x >= 0.5) {
  1140.     return lngamma_lanczos(x, result);
  1141.   }
  1142.   else if(x == 0.0) {
  1143.     DOMAIN_ERROR(result);
  1144.   }
  1145.   else if(fabs(x) < 0.02) {
  1146.     double sgn;
  1147.     return lngamma_sgn_0(x, result, &sgn);
  1148.   }
  1149.   else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) {
  1150.     /* Try to extract a fractional
  1151.      * part from x.
  1152.      */
  1153.     double z  = 1.0 - x;
  1154.     double s  = sin(M_PI*z);
  1155.     double as = fabs(s);
  1156.     if(s == 0.0) {
  1157.       DOMAIN_ERROR(result);
  1158.     }
  1159.     else if(as < M_PI*0.015) {
  1160.       /* x is near a negative integer, -N */
  1161.       if(x < INT_MIN + 2.0) {
  1162.         result->val = 0.0;
  1163.         result->err = 0.0;
  1164.         GSL_ERROR ("error", GSL_EROUND);
  1165.       }
  1166.       else {
  1167.         int N = -(int)(x - 0.5);
  1168.         double eps = x + N;
  1169.         double sgn;
  1170.         return lngamma_sgn_sing(N, eps, result, &sgn);
  1171.       }
  1172.     }
  1173.     else {
  1174.       gsl_sf_result lg_z;
  1175.       lngamma_lanczos(z, &lg_z);
  1176.       result->val = M_LNPI - (log(as) + lg_z.val);
  1177.       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg_z.err;
  1178.       return GSL_SUCCESS;
  1179.     }
  1180.   }
  1181.   else {
  1182.     /* |x| was too large to extract any fractional part */
  1183.     result->val = 0.0;
  1184.     result->err = 0.0;
  1185.     GSL_ERROR ("error", GSL_EROUND);
  1186.   }
  1187. }
  1188.  
  1189.  
  1190. int gsl_sf_lngamma_sgn_e(double x, gsl_sf_result * result_lg, double * sgn)
  1191. {
  1192.   if(fabs(x - 1.0) < 0.01) {
  1193.     *sgn = 1.0;
  1194.     return lngamma_1_pade(x-1.0, result_lg);
  1195.   }
  1196.   else if(fabs(x - 2.0) < 0.01) {
  1197.     *sgn = 1.0;
  1198.     return lngamma_2_pade(x-2.0, result_lg);
  1199.   }
  1200.   else if(x >= 0.5) {
  1201.     *sgn = 1.0;
  1202.     return lngamma_lanczos(x, result_lg);
  1203.   }
  1204.   else if(x == 0.0) {
  1205.     *sgn = 0.0;
  1206.     DOMAIN_ERROR(result_lg);
  1207.   }
  1208.   else if(fabs(x) < 0.02) {
  1209.     return lngamma_sgn_0(x, result_lg, sgn);
  1210.   }
  1211.   else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) {
  1212.     double z = 1.0 - x;
  1213.     double s = sin(M_PI*x);
  1214.     double as = fabs(s);
  1215.     if(s == 0.0) {
  1216.       *sgn = 0.0;
  1217.       DOMAIN_ERROR(result_lg);
  1218.     }
  1219.     else if(as < M_PI*0.015) {
  1220.       /* x is near a negative integer, -N */
  1221.       if(x < INT_MIN + 2.0) {
  1222.         result_lg->val = 0.0;
  1223.         result_lg->err = 0.0;
  1224.     *sgn = 0.0;
  1225.     GSL_ERROR ("error", GSL_EROUND);
  1226.       }
  1227.       else {
  1228.         int N = -(int)(x - 0.5);
  1229.     double eps = x + N;
  1230.         return lngamma_sgn_sing(N, eps, result_lg, sgn);
  1231.       }
  1232.     }
  1233.     else {
  1234.       gsl_sf_result lg_z;
  1235.       lngamma_lanczos(z, &lg_z);
  1236.       *sgn = (s > 0.0 ? 1.0 : -1.0);
  1237.       result_lg->val = M_LNPI - (log(as) + lg_z.val);
  1238.       result_lg->err = 2.0 * GSL_DBL_EPSILON * fabs(result_lg->val) + lg_z.err;
  1239.       return GSL_SUCCESS;
  1240.     }
  1241.   }
  1242.   else {
  1243.     /* |x| was too large to extract any fractional part */
  1244.     result_lg->val = 0.0;
  1245.     result_lg->err = 0.0;
  1246.     *sgn = 0.0;
  1247.     GSL_ERROR ("error", GSL_EROUND);
  1248.   }
  1249. }
  1250.  
  1251.  
  1252. int
  1253. gsl_sf_gamma_e(const double x, gsl_sf_result * result)
  1254. {
  1255.   if(x < 0.5) {
  1256.     int rint_x = (int)floor(x+0.5);
  1257.     double f_x = x - rint_x;
  1258.     double sgn = ( GSL_IS_EVEN(rint_x) ? 1.0 : -1.0 );
  1259.     double sin_term = sgn * sin(M_PI * f_x) / M_PI;
  1260.  
  1261.     if(sin_term == 0.0) {
  1262.       DOMAIN_ERROR(result);
  1263.     }
  1264.     else if(x > -169.0) {
  1265.       gsl_sf_result g;
  1266.       gamma_xgthalf(1.0-x, &g);
  1267.       if(fabs(sin_term) * g.val * GSL_DBL_MIN < 1.0) {
  1268.         result->val  = 1.0/(sin_term * g.val);
  1269.     result->err  = fabs(g.err/g.val) * fabs(result->val);
  1270.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1271.         return GSL_SUCCESS;
  1272.       }
  1273.       else {
  1274.         UNDERFLOW_ERROR(result);
  1275.       }
  1276.     }
  1277.     else {
  1278.       /* It is hard to control it here.
  1279.        * We can only exponentiate the
  1280.        * logarithm and eat the loss of
  1281.        * precision.
  1282.        */
  1283.       gsl_sf_result lng;
  1284.       double sgn;
  1285.       int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn);
  1286.       int stat_e   = gsl_sf_exp_mult_err_e(lng.val, lng.err, sgn, 0.0, result);
  1287.       return GSL_ERROR_SELECT_2(stat_e, stat_lng);
  1288.     }
  1289.   }
  1290.   else {
  1291.     return gamma_xgthalf(x, result);
  1292.   }
  1293. }
  1294.  
  1295.  
  1296. int
  1297. gsl_sf_gammastar_e(const double x, gsl_sf_result * result)
  1298. {
  1299.   /* CHECK_POINTER(result) */
  1300.  
  1301.   if(x <= 0.0) {
  1302.     DOMAIN_ERROR(result);
  1303.   }
  1304.   else if(x < 0.5) {
  1305.     gsl_sf_result lg;
  1306.     const int stat_lg = gsl_sf_lngamma_e(x, &lg);
  1307.     const double lx = log(x);
  1308.     const double c  = 0.5*(M_LN2+M_LNPI);
  1309.     const double lnr_val = lg.val - (x-0.5)*lx + x - c;
  1310.     const double lnr_err = lg.err + 2.0 * GSL_DBL_EPSILON *((x+0.5)*fabs(lx) + c);
  1311.     const int stat_e  = gsl_sf_exp_err_e(lnr_val, lnr_err, result);
  1312.     return GSL_ERROR_SELECT_2(stat_lg, stat_e);
  1313.   }
  1314.   else if(x < 2.0) {
  1315.     const double t = 4.0/3.0*(x-0.5) - 1.0;
  1316.     return cheb_eval_e(&gstar_a_cs, t, result);
  1317.   }
  1318.   else if(x < 10.0) {
  1319.     const double t = 0.25*(x-2.0) - 1.0;
  1320.     gsl_sf_result c;
  1321.     cheb_eval_e(&gstar_b_cs, t, &c);
  1322.     result->val  = c.val/(x*x) + 1.0 + 1.0/(12.0*x);
  1323.     result->err  = c.err/(x*x);
  1324.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1325.     return GSL_SUCCESS;
  1326.   }
  1327.   else if(x < 1.0/GSL_ROOT4_DBL_EPSILON) {
  1328.     return gammastar_ser(x, result);
  1329.   }
  1330.   else if(x < 1.0/GSL_DBL_EPSILON) {
  1331.     /* Use Stirling formula for Gamma(x).
  1332.      */
  1333.     const double xi = 1.0/x;
  1334.     result->val = 1.0 + xi/12.0*(1.0 + xi/24.0*(1.0 - xi*(139.0/180.0 + 571.0/8640.0*xi)));
  1335.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1336.     return GSL_SUCCESS;
  1337.   }
  1338.   else {
  1339.     result->val = 1.0;
  1340.     result->err = 1.0/x;
  1341.     return GSL_SUCCESS;
  1342.   }
  1343. }
  1344.  
  1345.  
  1346. int
  1347. gsl_sf_gammainv_e(const double x, gsl_sf_result * result)
  1348. {
  1349.   /* CHECK_POINTER(result) */
  1350.  
  1351.   if(x < 0.5) {
  1352.     gsl_sf_result lng;
  1353.     double sgn;
  1354.     int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn);
  1355.     if(stat_lng == GSL_EDOM) {
  1356.       result->val = 0.0;
  1357.       result->err = 0.0;
  1358.       return GSL_SUCCESS;
  1359.     }
  1360.     else if(stat_lng != GSL_SUCCESS) {
  1361.       result->val = 0.0;
  1362.       result->err = 0.0;
  1363.       return stat_lng;
  1364.     }
  1365.     else {
  1366.       return gsl_sf_exp_mult_err_e(-lng.val, lng.err, sgn, 0.0, result);
  1367.     }
  1368.   }
  1369.   else {
  1370.     gsl_sf_result g;
  1371.     int stat_g = gamma_xgthalf(x, &g);
  1372.     if(stat_g == GSL_EOVRFLW) {
  1373.       UNDERFLOW_ERROR(result);
  1374.     }
  1375.     else {
  1376.       result->val  = 1.0/g.val;
  1377.       result->err  = fabs(g.err/g.val) * fabs(result->val);
  1378.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1379.       CHECK_UNDERFLOW(result);
  1380.       return GSL_SUCCESS;
  1381.     }
  1382.   }
  1383. }
  1384.  
  1385.  
  1386. int
  1387. gsl_sf_lngamma_complex_e(double zr, double zi, gsl_sf_result * lnr, gsl_sf_result * arg)
  1388. {
  1389.   if(zr <= 0.5) {
  1390.     /* Transform to right half plane using reflection;
  1391.      * in fact we do a little better by stopping at 1/2.
  1392.      */
  1393.     double x = 1.0-zr;
  1394.     double y = -zi;
  1395.     gsl_sf_result a, b;
  1396.     gsl_sf_result lnsin_r, lnsin_i;
  1397.  
  1398.     int stat_l = lngamma_lanczos_complex(x, y, &a, &b);
  1399.     int stat_s = gsl_sf_complex_logsin_e(M_PI*zr, M_PI*zi, &lnsin_r, &lnsin_i);
  1400.  
  1401.     if(stat_s == GSL_SUCCESS) {
  1402.       int stat_r;
  1403.       lnr->val = M_LNPI - lnsin_r.val - a.val;
  1404.       lnr->err = lnsin_r.err + a.err + 2.0 * GSL_DBL_EPSILON * fabs(lnr->val);
  1405.       arg->val = -lnsin_i.val - b.val;
  1406.       arg->err = lnsin_i.err + b.err + 2.0 * GSL_DBL_EPSILON * fabs(arg->val);
  1407.       stat_r = gsl_sf_angle_restrict_symm_e(&(arg->val));
  1408.       return GSL_ERROR_SELECT_2(stat_r, stat_l);
  1409.     }
  1410.     else {
  1411.       DOMAIN_ERROR_2(lnr,arg);
  1412.     }
  1413.   }
  1414.   else {
  1415.     /* otherwise plain vanilla Lanczos */
  1416.     return lngamma_lanczos_complex(zr, zi, lnr, arg);
  1417.   }
  1418. }
  1419.  
  1420.  
  1421. int gsl_sf_taylorcoeff_e(const int n, const double x, gsl_sf_result * result)
  1422. {
  1423.   /* CHECK_POINTER(result) */
  1424.  
  1425.   if(x < 0.0 || n < 0) {
  1426.     DOMAIN_ERROR(result);
  1427.   }
  1428.   else if(n == 0) {
  1429.     result->val = 1.0;
  1430.     result->err = 0.0;
  1431.     return GSL_SUCCESS;
  1432.   }
  1433.   else if(n == 1) {
  1434.     result->val = x;
  1435.     result->err = 0.0;
  1436.     return GSL_SUCCESS;
  1437.   }
  1438.   else if(x == 0.0) {
  1439.     result->val = 0.0;
  1440.     result->err = 0.0;
  1441.     return GSL_SUCCESS;
  1442.   }
  1443.   else {
  1444.     const double log2pi = M_LNPI + M_LN2;
  1445.     const double ln_test = n*(log(x)+1.0) + 1.0 - (n+0.5)*log(n+1.0) + 0.5*log2pi;
  1446.  
  1447.     if(ln_test < GSL_LOG_DBL_MIN+1.0) {
  1448.       UNDERFLOW_ERROR(result);
  1449.     }
  1450.     else if(ln_test > GSL_LOG_DBL_MAX-1.0) {
  1451.       OVERFLOW_ERROR(result);
  1452.     }
  1453.     else {
  1454.       double product = 1.0;
  1455.       int k;
  1456.       for(k=1; k<=n; k++) {
  1457.         product *= (x/k);
  1458.       }
  1459.       result->val = product;
  1460.       result->err = n * GSL_DBL_EPSILON * product;
  1461.       CHECK_UNDERFLOW(result);
  1462.       return GSL_SUCCESS;
  1463.     }    
  1464.   }
  1465. }
  1466.  
  1467.  
  1468. int gsl_sf_fact_e(const unsigned int n, gsl_sf_result * result)
  1469. {
  1470.   /* CHECK_POINTER(result) */
  1471.  
  1472.   if(n < 18) {
  1473.     result->val = fact_table[n].f;
  1474.     result->err = 0.0;
  1475.     return GSL_SUCCESS;
  1476.   }
  1477.   else if(n <= FACT_TABLE_MAX){
  1478.     result->val = fact_table[n].f;
  1479.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1480.     return GSL_SUCCESS;
  1481.   }
  1482.   else {
  1483.     OVERFLOW_ERROR(result);
  1484.   }
  1485. }
  1486.  
  1487.  
  1488. int gsl_sf_doublefact_e(const unsigned int n, gsl_sf_result * result)
  1489. {
  1490.   /* CHECK_POINTER(result) */
  1491.  
  1492.   if(n < 26) {
  1493.     result->val = doub_fact_table[n].f;
  1494.     result->err = 0.0;
  1495.     return GSL_SUCCESS;
  1496.   }
  1497.   else if(n <= DOUB_FACT_TABLE_MAX){
  1498.     result->val = doub_fact_table[n].f;
  1499.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1500.     return GSL_SUCCESS;
  1501.   }
  1502.   else {
  1503.     OVERFLOW_ERROR(result);
  1504.   }
  1505. }
  1506.  
  1507.  
  1508. int gsl_sf_lnfact_e(const unsigned int n, gsl_sf_result * result)
  1509. {
  1510.   /* CHECK_POINTER(result) */
  1511.  
  1512.   if(n <= FACT_TABLE_MAX){
  1513.     result->val = log(fact_table[n].f);
  1514.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1515.     return GSL_SUCCESS;
  1516.   }
  1517.   else {
  1518.     gsl_sf_lngamma_e(n+1.0, result);
  1519.     return GSL_SUCCESS;
  1520.   }
  1521. }
  1522.  
  1523.  
  1524. int gsl_sf_lndoublefact_e(const unsigned int n, gsl_sf_result * result)
  1525. {
  1526.   /* CHECK_POINTER(result) */
  1527.  
  1528.   if(n <= DOUB_FACT_TABLE_MAX){
  1529.     result->val = log(doub_fact_table[n].f);
  1530.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1531.     return GSL_SUCCESS;
  1532.   }
  1533.   else if(GSL_IS_ODD(n)) {
  1534.     gsl_sf_result lg;
  1535.     gsl_sf_lngamma_e(0.5*(n+2.0), &lg);
  1536.     result->val = 0.5*(n+1.0) * M_LN2 - 0.5*M_LNPI + lg.val;
  1537.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err;
  1538.     return GSL_SUCCESS;
  1539.   }
  1540.   else {
  1541.     gsl_sf_result lg;
  1542.     gsl_sf_lngamma_e(0.5*n+1.0, &lg);
  1543.     result->val = 0.5*n*M_LN2 + lg.val;
  1544.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err;
  1545.     return GSL_SUCCESS;
  1546.   }
  1547. }
  1548.  
  1549.  
  1550. int gsl_sf_lnchoose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
  1551. {
  1552.   /* CHECK_POINTER(result) */
  1553.  
  1554.   if(m > n) {
  1555.     DOMAIN_ERROR(result);
  1556.   }
  1557.   else if(m == n || m == 0) {
  1558.     result->val = 0.0;
  1559.     result->err = 0.0;
  1560.     return GSL_SUCCESS;
  1561.   }
  1562.   else {
  1563.     gsl_sf_result nf;
  1564.     gsl_sf_result mf;
  1565.     gsl_sf_result nmmf;
  1566.     if(m*2 > n) m = n-m;
  1567.     gsl_sf_lnfact_e(n, &nf);
  1568.     gsl_sf_lnfact_e(m, &mf);
  1569.     gsl_sf_lnfact_e(n-m, &nmmf);
  1570.     result->val  = nf.val - mf.val - nmmf.val;
  1571.     result->err  = nf.err + mf.err + nmmf.err;
  1572.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1573.     return GSL_SUCCESS;
  1574.   }
  1575. }
  1576.  
  1577.  
  1578. int gsl_sf_choose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
  1579. {
  1580.   if(m > n) {
  1581.     DOMAIN_ERROR(result);
  1582.   }
  1583.   else if(m == n || m == 0) {
  1584.     result->val = 1.0;
  1585.     result->err = 0.0;
  1586.     return GSL_SUCCESS;
  1587.   }
  1588.   else {
  1589.     double prod = 1.0;
  1590.     int k;
  1591.     for(k=n; k>=m+1; k--) {
  1592.       double tk = (double)k / (double)(k-m);
  1593.       if(tk > GSL_DBL_MAX/prod) {
  1594.         OVERFLOW_ERROR(result);
  1595.       }
  1596.       prod *= tk;
  1597.     }
  1598.     result->val = prod;
  1599.     result->err = 2.0 * GSL_DBL_EPSILON * prod * fabs(n-m);
  1600.     return GSL_SUCCESS;
  1601.   }
  1602. }
  1603.  
  1604.  
  1605. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  1606.  
  1607. #include "eval.h"
  1608.  
  1609. double gsl_sf_fact(const unsigned int n)
  1610. {
  1611.   EVAL_RESULT(gsl_sf_fact_e(n, &result));
  1612. }
  1613.  
  1614. double gsl_sf_lnfact(const unsigned int n)
  1615. {
  1616.   EVAL_RESULT(gsl_sf_lnfact_e(n, &result));
  1617. }
  1618.  
  1619. double gsl_sf_doublefact(const unsigned int n)
  1620. {
  1621.   EVAL_RESULT(gsl_sf_doublefact_e(n, &result));
  1622. }
  1623.  
  1624. double gsl_sf_lndoublefact(const unsigned int n)
  1625. {
  1626.   EVAL_RESULT(gsl_sf_lndoublefact_e(n, &result));
  1627. }
  1628.  
  1629. double gsl_sf_lngamma(const double x)
  1630. {
  1631.   EVAL_RESULT(gsl_sf_lngamma_e(x, &result));
  1632. }
  1633.  
  1634. double gsl_sf_gamma(const double x)
  1635. {
  1636.   EVAL_RESULT(gsl_sf_gamma_e(x, &result));
  1637. }
  1638.  
  1639. double gsl_sf_gammastar(const double x)
  1640. {
  1641.   EVAL_RESULT(gsl_sf_gammastar_e(x, &result));
  1642. }
  1643.  
  1644. double gsl_sf_gammainv(const double x)
  1645. {
  1646.   EVAL_RESULT(gsl_sf_gammainv_e(x, &result));
  1647. }
  1648.  
  1649. double gsl_sf_taylorcoeff(const int n, const double x)
  1650. {
  1651.   EVAL_RESULT(gsl_sf_taylorcoeff_e(n, x, &result));
  1652. }
  1653.  
  1654. double gsl_sf_choose(unsigned int n, unsigned int m)
  1655. {
  1656.   EVAL_RESULT(gsl_sf_choose_e(n, m, &result));
  1657. }
  1658.  
  1659. double gsl_sf_lnchoose(unsigned int n, unsigned int m)
  1660. {
  1661.   EVAL_RESULT(gsl_sf_lnchoose_e(n, m, &result));
  1662. }
  1663.